##########################################
### LANE AND SCHOENHERR CITE AND SWAY? ###
##########################################

# Purpose:
# 	- Replicate Table 4
# 	- Prepare data for Figures 2 and 3, produced in Analysis2Plots.R

# Data Sources:
#	- Analysis2Data.csv

# Run on:
#	- R 4.4.0 ("Puppy Cup") on macOS Sequoia 15.1.1.

# Associated codebook:
#	- Analysis2Codebook.pdf

###########################################
###########################################
###########################################

library(foreign)
library(arm)
library(readxl)
library(lattice)
library(tidyverse)
library(stargazer)
library(faraway)
library(lme4)

######################################
### STEP 0: SET WORKDING DIRECTORY ###
######################################

setwd()

########################################
### STEP 1: PULL IN THE MODEL 2 DATA ###
########################################

model2data <- read.csv("Analysis2Data.csv")

############################################
### STEP 2: PRE-PROCESS SOME OF THE DATA ###
############################################

model2data$logPetBriefExperience <- log(model2data$petBriefExperience + 1)
model2data$logRespBriefExperience <- log(model2data$respBriefExperience + 1)

model2data$logPetNumCites <- log(model2data$petNumCites + 1)
model2data$logRespNumCites <- log(model2data$respNumCites + 1)

# check all the variables to make sure no NAs before running:
table(is.na(model2data$voteWithPet))
table(is.na(model2data$petUnnamedCiteCount))
table(is.na(model2data$petMajCallOutCount))
table(is.na(model2data$petNotMajOpinCount))
table(is.na(model2data$respUnnamedCiteCount))
table(is.na(model2data$respMajCallOutCount))
table(is.na(model2data$respNotMajOpinCount))
table(is.na(model2data$ideoAlign))
table(is.na(model2data$logPetBriefExperience))
table(is.na(model2data$logRespBriefExperience))
table(is.na(model2data$sgParty))
table(is.na(model2data$helpOSG))
table(is.na(model2data$amiciNet))
table(is.na(model2data$petFormerClerk))
table(is.na(model2data$respFormerClerk))
table(is.na(model2data$netStatus))
table(is.na(model2data$lcDisagreement))
model2data <- model2data %>% filter(!is.na(lcDisagreement))
table(is.na(model2data$CSI))
table(is.na(model2data$pastExpertise))
table(is.na(model2data$petNumCites))
table(is.na(model2data$respNumCites))
table(is.na(model2data$petQuestionsOA))
table(is.na(model2data$respQuestionsOA))
table(is.na(model2data$justice))
table(is.na(model2data$issueArea))
table(is.na(model2data$term))

# convert the three levels to factors
model2data$justice <- as.factor(model2data$justice)
model2data$issueArea <- as.factor(model2data$issueArea)
model2data$term <- as.factor(model2data$term)

#####################
### STEP 3: MODEL ###
#####################

### TABLE 4 ###

model2data$petUnnamedCiteLog <- log(model2data$petUnnamedCiteCount + 1)
model2data$petMajCallOutLog <- log(model2data$petMajCallOutCount + 1)
model2data$petNotMajOpinLog <- log(model2data$petNotMajOpinCount + 1)
model2data$respUnnamedCiteLog <- log(model2data$respUnnamedCiteCount + 1)
model2data$respMajCallOutLog <- log(model2data$respMajCallOutCount + 1)
model2data$respNotMajOpinLog <- log(model2data$respNotMajOpinCount + 1)
model2data$pastExpertiseLog <- log(model2data$pastExpertise + 1)
model2data$oaQuestDiff <- model2data$petQuestionsOA - model2data$respQuestionsOA
model2data$petExperienceAdvantage <- model2data$petBriefExperience - model2data$respBriefExperience

model2logCites <- glmer(voteWithPet ~ petUnnamedCiteLog
							+ petMajCallOutLog
							+ petNotMajOpinLog
							+ respUnnamedCiteLog
							+ respMajCallOutLog
							+ respNotMajOpinLog
							+ ideoAlign
							+ pastExpertise
							+ logPetNumCites
							+ logRespNumCites
							+ petExperienceAdvantage
							+ sgParty
							+ lcDisagreement
							+ amiciNet
							+ helpOSG
							+ netStatus
							+ oaQuestDiff
							+ petUnnamedCiteLog * ideoAlign
							+ petMajCallOutLog * ideoAlign
							+ petNotMajOpinLog * ideoAlign
							+ respUnnamedCiteLog * ideoAlign
							+ respMajCallOutLog * ideoAlign
							+ respNotMajOpinLog * ideoAlign
							+ (1 | issueArea)
							+ (1 | term)
							+ (1 | justice),
							data = model2data,
							family = binomial(link = logit),
							glmerControl(optimizer = "Nelder_Mead"))

summary(model2logCites)
display(model2logCites)

####################################
### STEP 3B: WRITE OUT THE TABLE ###
####################################

stargazer(model2logCites, align = TRUE, omit.stat=c("LL", "ser", "f"))

###################################
### STEP 3C: KEEP FOR REFERENCE ###
###################################

# summary statistics
summary(model2data$petUnnamedCiteLog)
summary(model2data$petMajCallOutLog)
summary(model2data$petNotMajOpinLog)
summary(model2data$respUnnamedCiteLog)
summary(model2data$respMajCallOutLog)
summary(model2data$respNotMajOpinLog)

################################################
### STEP 4: PREP FOR PREDICTED PROBABILITIES ###
################################################

### create the datasets for predProbs for the different intercepts
model2data <- model2data %>% arrange(term, justice, issueArea)
model2data$count <- 1

# term
term <- model2data %>% pivot_wider(names_from = term, values_from = count, values_fill = 0)
term <- term[134:168]

# issueArea
issueArea <- model2data %>% pivot_wider(names_from = issueArea, values_from = count, values_fill = 0)
issueArea <- issueArea[134:145]

# justice
justice <- model2data %>% pivot_wider(names_from = justice, values_from = count, values_fill = 0)
justice <- justice[134:154]

#######################################
### STEP 5: PREDICTED PROBABILITIES ###
#######################################

# values needed: (0.05 = 0, 0.95 = 11)
# remember the converstion is ln(x+1) to avoid indefinite answers
#	- numbers below are x, not x+1, but you need x+1 to calculate
# 0 - 0
# 1 - 0.693147
# 2 - 1.098612
# 3 - 1.386294
# 4 - 1.609438
# 5 - 1.791959
# 6 - 1.945910
# 7 - 2.079442
# 8 - 2.197225
# 9 - 2.302585
# 10 - 2.397895
# 11 - 2.484907

### petUnnamedCiteLog x ideoAlign
# low = quantile(model2data$ideoAlign, 0.1) = -2.909
# medium = quantile(model2data$ideoAlign, 0.5) = 0.077
# high = quantile(model2data$ideoAlign, 0.9) = 2.918
# this is Thomas/Ginsburg area
#	- so align is Ginsburg + Planned Parenthood or Thomas and the NRA
#	- Not aligned is Ginsburg + NRA, Thomas + Planned Parenthood

##################################
### STEP 5A: PET PASSIVE CITES ###
##################################

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model2logCites, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

#petPassiveSeq <- seq(from = quantile(model2data$petUnnamedCiteLog, 0.05),
#					to = quantile(model2data$petUnnamedCiteLog, 0.95),
#					by = 0.1)

petPassiveSeq <- c(0, 0.693147, 1.098612, 1.386294, 1.609438, 1.791959, 1.945910, 2.079442, 2.197225, 2.302585, 2.397895, 2.484907)
					
pMeanPetPassiveLowCong <- as.matrix(NA, nrow = length(petPassiveSeq))
pLoPetPassiveLowCong <- as.matrix(NA, nrow = length(petPassiveSeq))
pHiPetPassiveLowCong <- as.matrix(NA, nrow = length(petPassiveSeq))

pMeanPetPassiveMedCong <- as.matrix(NA, nrow = length(petPassiveSeq))
pLoPetPassiveMedCong <- as.matrix(NA, nrow = length(petPassiveSeq))
pHiPetPassiveMedCong <- as.matrix(NA, nrow = length(petPassiveSeq))

pMeanPetPassiveHighCong <- as.matrix(NA, nrow = length(petPassiveSeq))
pLoPetPassiveHighCong <- as.matrix(NA, nrow = length(petPassiveSeq))
pHiPetPassiveHighCong <- as.matrix(NA, nrow = length(petPassiveSeq))

pMeanPetPassiveDiffLowHigh <- as.matrix(NA, nrow = length(petPassiveSeq))
pLoPetPassiveDiffLowHigh <- as.matrix(NA, nrow = length(petPassiveSeq))
pHiPetPassiveDiffLowHigh <- as.matrix(NA, nrow = length(petPassiveSeq))

pMeanPetPassiveDiffLowMed <- as.matrix(NA, nrow = length(petPassiveSeq))
pLoPetPassiveDiffLowMed <- as.matrix(NA, nrow = length(petPassiveSeq))
pHiPetPassiveDiffLowMed <- as.matrix(NA, nrow = length(petPassiveSeq))

pMeanPetPassiveDiffMedHigh <- as.matrix(NA, nrow = length(petPassiveSeq))
pLoPetPassiveDiffMedHigh <- as.matrix(NA, nrow = length(petPassiveSeq))
pHiPetPassiveDiffMedHigh <- as.matrix(NA, nrow = length(petPassiveSeq))

for(i in 1:length(petPassiveSeq)){
	predPetPassiveLowCong <- cbind(1,
							petPassiveSeq[i],
							model2data$petMajCallOutLog,
							model2data$petNotMajOpinLog,
							model2data$respUnnamedCiteLog,
							model2data$respMajCallOutLog,
							model2data$respNotMajOpinLog,
							-2.909,
							model2data$pastExpertise,
							model2data$logPetNumCites,
							model2data$logRespNumCites,
							model2data$petExperienceAdvantage,
							model2data$sgParty,
							model2data$lcDisagreement,
							model2data$amiciNet,
							model2data$helpOSG,
							model2data$netStatus,
							model2data$oaQuestDiff,
							petPassiveSeq[i] * -2.909,
							model2data$petMajCallOutLog * -2.909,
							model2data$petNotMajOpinLog * -2.909,
							model2data$respUnnamedCiteLog * -2.909,
							model2data$respMajCallOutLog * -2.909,
							model2data$respNotMajOpinLog * -2.909,
							issueArea,
							term,
							justice
							)
	predPetPassiveLowCong <- as.matrix(predPetPassiveLowCong)
	
	predPetPassiveMedCong <- cbind(1,
							petPassiveSeq[i],
							model2data$petMajCallOutLog,
							model2data$petNotMajOpinLog,
							model2data$respUnnamedCiteLog,
							model2data$respMajCallOutLog,
							model2data$respNotMajOpinLog,
							0.077,
							model2data$pastExpertise,
							model2data$logPetNumCites,
							model2data$logRespNumCites,	
							model2data$petExperienceAdvantage,
							model2data$sgParty,
							model2data$lcDisagreement,
							model2data$amiciNet,
							model2data$helpOSG,
							model2data$netStatus,						
							model2data$oaQuestDiff,
							petPassiveSeq[i] * 0.077,
							model2data$petMajCallOutLog * 0.077,
							model2data$petNotMajOpinLog * 0.077,
							model2data$respUnnamedCiteLog * 0.077,
							model2data$respMajCallOutLog * 0.077,
							model2data$respNotMajOpinLog * 0.077,
							issueArea,
							term,
							justice
							)
	predPetPassiveMedCong <- as.matrix(predPetPassiveMedCong)
	
	predPetPassiveHighCong <- cbind(1,
							petPassiveSeq[i],
							model2data$petMajCallOutLog,
							model2data$petNotMajOpinLog,
							model2data$respUnnamedCiteLog,
							model2data$respMajCallOutLog,
							model2data$respNotMajOpinLog,
							2.918,
							model2data$pastExpertise,
							model2data$logPetNumCites,
							model2data$logRespNumCites,	
							model2data$petExperienceAdvantage,
							model2data$sgParty,
							model2data$lcDisagreement,
							model2data$amiciNet,
							model2data$helpOSG,
							model2data$netStatus,						
							model2data$oaQuestDiff,
							petPassiveSeq[i] * 2.918,
							model2data$petMajCallOutLog * 2.918,
							model2data$petNotMajOpinLog * 2.918,
							model2data$respUnnamedCiteLog * 2.918,
							model2data$respMajCallOutLog * 2.918,
							model2data$respNotMajOpinLog * 2.918,
							issueArea,
							term,
							justice
							)
	predPetPassiveHighCong <- as.matrix(predPetPassiveHighCong)
	
	ppPetPassiveLowCong <- as.matrix(NA, row = n.draws)
	ppPetPassiveMedCong <- as.matrix(NA, row = n.draws)
	ppPetPassiveHighCong <- as.matrix(NA, row = n.draws)
	ppPetPassiveDiffLowHigh <- as.matrix(NA, row = n.draws)
	ppPetPassiveDiffLowMed <- as.matrix(NA, row = n.draws)
	ppPetPassiveDiffMedHigh <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbPetPassiveLowCong <- ilogit(predPetPassiveLowCong %*% t(sim.coef[j, ]))
		ppPetPassiveLowCong[j] <- mean(xbPetPassiveLowCong)
		xbPetPassiveMedCong <- ilogit(predPetPassiveMedCong %*% t(sim.coef[j, ]))
		ppPetPassiveMedCong[j] <- mean(xbPetPassiveMedCong)
		xbPetPassiveHighCong <- ilogit(predPetPassiveHighCong %*% t(sim.coef[j, ]))
		ppPetPassiveHighCong[j] <- mean(xbPetPassiveHighCong)
		ppPetPassiveDiffLowHigh[j] <- mean(xbPetPassiveHighCong - xbPetPassiveLowCong)
		ppPetPassiveDiffLowMed[j] <- mean(xbPetPassiveMedCong - xbPetPassiveLowCong)
		ppPetPassiveDiffMedHigh[j] <- mean(xbPetPassiveHighCong - xbPetPassiveMedCong)
		
	}
	
	pMeanPetPassiveLowCong[i] <- mean(ppPetPassiveLowCong)
	pLoPetPassiveLowCong[i] <- quantile(ppPetPassiveLowCong, 0.025, na.rm = TRUE)
	pHiPetPassiveLowCong[i] <- quantile(ppPetPassiveLowCong, 0.975, na.rm = TRUE)
	
	pMeanPetPassiveMedCong[i] <- mean(ppPetPassiveMedCong)
	pLoPetPassiveMedCong[i] <- quantile(ppPetPassiveMedCong, 0.025, na.rm = TRUE)
	pHiPetPassiveMedCong[i] <- quantile(ppPetPassiveMedCong, 0.975, na.rm = TRUE)
	
	pMeanPetPassiveHighCong[i] <- mean(ppPetPassiveHighCong)
	pLoPetPassiveHighCong[i] <- quantile(ppPetPassiveHighCong, 0.025, na.rm = TRUE)
	pHiPetPassiveHighCong[i] <- quantile(ppPetPassiveHighCong, 0.975, na.rm = TRUE)
	
	pMeanPetPassiveDiffLowHigh[i] <- mean(ppPetPassiveDiffLowHigh)
	pLoPetPassiveDiffLowHigh[i] <- quantile(ppPetPassiveDiffLowHigh, 0.025, na.rm = TRUE)
	pHiPetPassiveDiffLowHigh[i] <- quantile(ppPetPassiveDiffLowHigh, 0.975, na.rm = TRUE)
	
	pMeanPetPassiveDiffLowMed[i] <- mean(ppPetPassiveDiffLowMed)
	pLoPetPassiveDiffLowMed[i] <- quantile(ppPetPassiveDiffLowMed, 0.025, na.rm = TRUE)
	pHiPetPassiveDiffLowMed[i] <- quantile(ppPetPassiveDiffLowMed, 0.975, na.rm = TRUE)
	
	pMeanPetPassiveDiffMedHigh[i] <- mean(ppPetPassiveDiffMedHigh)
	pLoPetPassiveDiffMedHigh[i] <- quantile(ppPetPassiveDiffMedHigh, 0.025, na.rm = TRUE)
	pHiPetPassiveDiffMedHigh[i] <- quantile(ppPetPassiveDiffMedHigh, 0.975, na.rm = TRUE)
	
}

cbind(pLoPetPassiveDiffLowHigh, pMeanPetPassiveDiffLowHigh, pHiPetPassiveDiffLowHigh)
cbind(pLoPetPassiveDiffLowMed, pMeanPetPassiveDiffLowMed, pHiPetPassiveDiffLowMed)
cbind(pLoPetPassiveDiffMedHigh, pMeanPetPassiveDiffMedHigh, pHiPetPassiveDiffMedHigh)

counter <- seq(0, 11, by = 1)

petCiteXcong <- data.frame(counter, petPassiveSeq, pLoPetPassiveLowCong, pMeanPetPassiveLowCong, pHiPetPassiveLowCong, pLoPetPassiveHighCong, pMeanPetPassiveHighCong, pHiPetPassiveHighCong, pLoPetPassiveMedCong, pMeanPetPassiveMedCong, pHiPetPassiveMedCong)
colnames(petCiteXcong) <- c("count", "passiveCite", "lbLow", "estLow", "ubLow", "lbHigh", "estHigh", "ubHigh", "lbMed", "estMed", "ubMed")
petCiteXcong$pet <- 1

# just for reference - this is not the table in the paper
#petPassiveCiteTable <- ggplot(petCiteXcong, aes(x = passiveCite)) + 
#				geom_line(aes(y = estLow), color = "red") +
#				geom_line(aes(y = lbLow), linetype = "dashed", color = "red") +
#				geom_line(aes(y = ubLow), linetype = "dashed", color = "red") +
#				geom_line(aes(y = estHigh), color = "blue") +
#				geom_line(aes(y = lbHigh), linetype = "dashed", color = "blue") +
#				geom_line(aes(y = ubHigh), linetype = "dashed", color = "blue") +
#				geom_line(aes(y = estMed), color = "green") +
#				geom_line(aes(y = lbMed), linetype = "dashed", color = "green") +
#				geom_line(aes(y = ubMed), linetype = "dashed", color = "green") + 
#				theme_light()								
#petPassiveCiteTable

###################################
### STEP 5B: RESP PASSIVE CITES ###
###################################

### petNotMajOpinLog x ideoAlign
# low = quantile(model2data$ideoAlign, 0.1) = -2.909
# medium = quantile(model2data$ideoAlign, 0.5) = 0.077
# high = quantile(model2data$ideoAlign, 0.9) = 2.918
# sequence = ln(count + 1), translates to 0-11

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model2logCites, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

respPassiveSeq <- c(0, 0.693147, 1.098612, 1.386294, 1.609438, 1.791959, 1.945910, 2.079442, 2.197225, 2.302585, 2.397895, 2.484907)
					
pMeanRespPassiveLowCong <- as.matrix(NA, nrow = length(respPassiveSeq))
pLoRespPassiveLowCong <- as.matrix(NA, nrow = length(respPassiveSeq))
pHiRespPassiveLowCong <- as.matrix(NA, nrow = length(respPassiveSeq))

pMeanRespPassiveMedCong <- as.matrix(NA, nrow = length(respPassiveSeq))
pLoRespPassiveMedCong <- as.matrix(NA, nrow = length(respPassiveSeq))
pHiRespPassiveMedCong <- as.matrix(NA, nrow = length(respPassiveSeq))

pMeanRespPassiveHighCong <- as.matrix(NA, nrow = length(respPassiveSeq))
pLoRespPassiveHighCong <- as.matrix(NA, nrow = length(respPassiveSeq))
pHiRespPassiveHighCong <- as.matrix(NA, nrow = length(respPassiveSeq))

pMeanRespPassiveDiffLowHigh <- as.matrix(NA, nrow = length(respPassiveSeq))
pLoRespPassiveDiffLowHigh <- as.matrix(NA, nrow = length(respPassiveSeq))
pHiRespPassiveDiffLowHigh <- as.matrix(NA, nrow = length(respPassiveSeq))

pMeanRespPassiveDiffLowMed <- as.matrix(NA, nrow = length(respPassiveSeq))
pLoRespPassiveDiffLowMed <- as.matrix(NA, nrow = length(respPassiveSeq))
pHiRespPassiveDiffLowMed <- as.matrix(NA, nrow = length(respPassiveSeq))

pMeanRespPassiveDiffMedHigh <- as.matrix(NA, nrow = length(respPassiveSeq))
pLoRespPassiveDiffMedHigh <- as.matrix(NA, nrow = length(respPassiveSeq))
pHiRespPassiveDiffMedHigh <- as.matrix(NA, nrow = length(respPassiveSeq))

for(i in 1:length(respPassiveSeq)){
	predRespPassiveLowCong <- cbind(1,
							model2data$petUnnamedCiteLog,
							model2data$petMajCallOutLog,
							model2data$petNotMajOpinLog,
							respPassiveSeq[i],
							model2data$respMajCallOutLog,
							model2data$respNotMajOpinLog,
							-2.909,
							model2data$pastExpertise,
							model2data$logPetNumCites,
							model2data$logRespNumCites,
							model2data$petExperienceAdvantage,
							model2data$sgParty,
							model2data$lcDisagreement,
							model2data$amiciNet,
							model2data$helpOSG,
							model2data$netStatus,
							model2data$oaQuestDiff,
							model2data$petUnnamedCiteLog * -2.909,
							model2data$petMajCallOutLog * -2.909,
							model2data$petNotMajOpinLog * -2.909,
							respPassiveSeq[i] * -2.909,
							model2data$respMajCallOutLog * -2.909,
							model2data$respNotMajOpinLog * -2.909,
							issueArea,
							term,
							justice
							)
	predRespPassiveLowCong <- as.matrix(predRespPassiveLowCong)
	
	predRespPassiveMedCong <- cbind(1,
							model2data$petUnnamedCiteLog,
							model2data$petMajCallOutLog,
							model2data$petNotMajOpinLog,
							respPassiveSeq[i],
							model2data$respMajCallOutLog,
							model2data$respNotMajOpinLog,
							0.077,
							model2data$pastExpertise,
							model2data$logPetNumCites,
							model2data$logRespNumCites,	
							model2data$petExperienceAdvantage,
							model2data$sgParty,
							model2data$lcDisagreement,
							model2data$amiciNet,
							model2data$helpOSG,
							model2data$netStatus,						
							model2data$oaQuestDiff,
							model2data$petUnnamedCiteLog * 0.077,
							model2data$petMajCallOutLog * 0.077,
							model2data$petNotMajOpinLog * 0.077,
							respPassiveSeq[i] * 0.077,
							model2data$respMajCallOutLog * 0.077,
							model2data$respNotMajOpinLog * 0.077,
							issueArea,
							term,
							justice
							)
	predRespPassiveMedCong <- as.matrix(predRespPassiveMedCong)
	
	predRespPassiveHighCong <- cbind(1,
							model2data$petUnnamedCiteLog,
							model2data$petMajCallOutLog,
							model2data$petNotMajOpinLog,
							respPassiveSeq[i],
							model2data$respMajCallOutLog,
							model2data$respNotMajOpinLog,
							2.918,
							model2data$pastExpertise,
							model2data$logPetNumCites,
							model2data$logRespNumCites,	
							model2data$petExperienceAdvantage,
							model2data$sgParty,
							model2data$lcDisagreement,
							model2data$amiciNet,
							model2data$helpOSG,
							model2data$netStatus,						
							model2data$oaQuestDiff,
							model2data$petUnnamedCiteLog * 2.918,
							model2data$petMajCallOutLog * 2.918,
							model2data$petNotMajOpinLog * 2.918,
							respPassiveSeq[i] * 2.918,
							model2data$respMajCallOutLog * 2.918,
							model2data$respNotMajOpinLog * 2.918,
							issueArea,
							term,
							justice
							)
	predRespPassiveHighCong <- as.matrix(predRespPassiveHighCong)
	
	ppRespPassiveLowCong <- as.matrix(NA, row = n.draws)
	ppRespPassiveMedCong <- as.matrix(NA, row = n.draws)
	ppRespPassiveHighCong <- as.matrix(NA, row = n.draws)
	ppRespPassiveDiffLowHigh <- as.matrix(NA, row = n.draws)
	ppRespPassiveDiffLowMed <- as.matrix(NA, row = n.draws)
	ppRespPassiveDiffMedHigh <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbRespPassiveLowCong <- ilogit(predRespPassiveLowCong %*% t(sim.coef[j, ]))
		ppRespPassiveLowCong[j] <- mean(xbRespPassiveLowCong)
		xbRespPassiveMedCong <- ilogit(predRespPassiveMedCong %*% t(sim.coef[j, ]))
		ppRespPassiveMedCong[j] <- mean(xbRespPassiveMedCong)
		xbRespPassiveHighCong <- ilogit(predRespPassiveHighCong %*% t(sim.coef[j, ]))
		ppRespPassiveHighCong[j] <- mean(xbRespPassiveHighCong)
		ppRespPassiveDiffLowHigh[j] <- mean(xbRespPassiveHighCong - xbRespPassiveLowCong)
		ppRespPassiveDiffLowMed[j] <- mean(xbRespPassiveMedCong - xbRespPassiveLowCong)
		ppRespPassiveDiffMedHigh[j] <- mean(xbRespPassiveHighCong - xbRespPassiveMedCong)
		
	}
	
	pMeanRespPassiveLowCong[i] <- mean(ppRespPassiveLowCong)
	pLoRespPassiveLowCong[i] <- quantile(ppRespPassiveLowCong, 0.025, na.rm = TRUE)
	pHiRespPassiveLowCong[i] <- quantile(ppRespPassiveLowCong, 0.975, na.rm = TRUE)
	
	pMeanRespPassiveMedCong[i] <- mean(ppRespPassiveMedCong)
	pLoRespPassiveMedCong[i] <- quantile(ppRespPassiveMedCong, 0.025, na.rm = TRUE)
	pHiRespPassiveMedCong[i] <- quantile(ppRespPassiveMedCong, 0.975, na.rm = TRUE)
	
	pMeanRespPassiveHighCong[i] <- mean(ppRespPassiveHighCong)
	pLoRespPassiveHighCong[i] <- quantile(ppRespPassiveHighCong, 0.025, na.rm = TRUE)
	pHiRespPassiveHighCong[i] <- quantile(ppRespPassiveHighCong, 0.975, na.rm = TRUE)
	
	pMeanRespPassiveDiffLowHigh[i] <- mean(ppRespPassiveDiffLowHigh)
	pLoRespPassiveDiffLowHigh[i] <- quantile(ppRespPassiveDiffLowHigh, 0.025, na.rm = TRUE)
	pHiRespPassiveDiffLowHigh[i] <- quantile(ppRespPassiveDiffLowHigh, 0.975, na.rm = TRUE)
	
	pMeanRespPassiveDiffLowMed[i] <- mean(ppRespPassiveDiffLowMed)
	pLoRespPassiveDiffLowMed[i] <- quantile(ppRespPassiveDiffLowMed, 0.025, na.rm = TRUE)
	pHiRespPassiveDiffLowMed[i] <- quantile(ppRespPassiveDiffLowMed, 0.975, na.rm = TRUE)
	
	pMeanRespPassiveDiffMedHigh[i] <- mean(ppRespPassiveDiffMedHigh)
	pLoRespPassiveDiffMedHigh[i] <- quantile(ppRespPassiveDiffMedHigh, 0.025, na.rm = TRUE)
	pHiRespPassiveDiffMedHigh[i] <- quantile(ppRespPassiveDiffMedHigh, 0.975, na.rm = TRUE)
	
}

cbind(pLoRespPassiveDiffLowHigh, pMeanRespPassiveDiffLowHigh, pHiRespPassiveDiffLowHigh)
cbind(pLoRespPassiveDiffLowMed, pMeanRespPassiveDiffLowMed, pHiRespPassiveDiffLowMed)
cbind(pLoRespPassiveDiffMedHigh, pMeanRespPassiveDiffMedHigh, pHiRespPassiveDiffMedHigh)

counter <- seq(0, 11, by = 1)

respCiteXcong <- data.frame(counter, respPassiveSeq, pLoRespPassiveLowCong, pMeanRespPassiveLowCong, pHiRespPassiveLowCong, pLoRespPassiveHighCong, pMeanRespPassiveHighCong, pHiRespPassiveHighCong, pLoRespPassiveMedCong, pMeanRespPassiveMedCong, pHiRespPassiveMedCong)
colnames(respCiteXcong) <- c("count", "passiveCite", "lbLow", "estLow", "ubLow", "lbHigh", "estHigh", "ubHigh", "lbMed", "estMed", "ubMed")
respCiteXcong$pet <- 0

# just for reference - this is not the table in the paper
respPassiveCiteTable <- ggplot(respCiteXcong, aes(x = passiveCite)) + 
				geom_line(aes(y = estLow), color = "red") +
				geom_line(aes(y = lbLow), linetype = "dashed", color = "red") +
				geom_line(aes(y = ubLow), linetype = "dashed", color = "red") +
				geom_line(aes(y = estHigh), color = "blue") +
				geom_line(aes(y = lbHigh), linetype = "dashed", color = "blue") +
				geom_line(aes(y = ubHigh), linetype = "dashed", color = "blue") +
				geom_line(aes(y = estMed), color = "green") +
				geom_line(aes(y = lbMed), linetype = "dashed", color = "green") +
				geom_line(aes(y = ubMed), linetype = "dashed", color = "green") + 
				theme_light()								
respPassiveCiteTable

#####################################
### STEP 5C: WRITE OUT FOR GRAPHS ###
#####################################

citesXcong <- rbind(petCiteXcong, respCiteXcong)

write.csv(citesXcong, "Analysis2PassivePredProbs.csv")

##################################
### STEP 6: SEPARATE CITATIONS ###
##################################

### petNotMajOpinLog x ideoAlign
# low = quantile(model2data$ideoAlign, 0.1) = -2.909
# medium = quantile(model2data$ideoAlign, 0.5) = 0.077
# high = quantile(model2data$ideoAlign, 0.9) = 2.918
# 0 - 0
# 1 - 0.693147
# 2 - 1.098612

########################
### STEP 6A: PET SEP ###
########################

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model2logCites, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

petSeparateSeq <- c(0, 0.693147, 1.098612)
					
pMeanPetSeparateLowCong <- as.matrix(NA, nrow = length(petSeparateSeq))
pLoPetSeparateLowCong <- as.matrix(NA, nrow = length(petSeparateSeq))
pHiPetSeparateLowCong <- as.matrix(NA, nrow = length(petSeparateSeq))

pMeanPetSeparateMedCong <- as.matrix(NA, nrow = length(petSeparateSeq))
pLoPetSeparateMedCong <- as.matrix(NA, nrow = length(petSeparateSeq))
pHiPetSeparateMedCong <- as.matrix(NA, nrow = length(petSeparateSeq))

pMeanPetSeparateHighCong <- as.matrix(NA, nrow = length(petSeparateSeq))
pLoPetSeparateHighCong <- as.matrix(NA, nrow = length(petSeparateSeq))
pHiPetSeparateHighCong <- as.matrix(NA, nrow = length(petSeparateSeq))

pMeanPetSeparateDiffLowHigh <- as.matrix(NA, nrow = length(petSeparateSeq))
pLoPetSeparateDiffLowHigh <- as.matrix(NA, nrow = length(petSeparateSeq))
pHiPetSeparateDiffLowHigh <- as.matrix(NA, nrow = length(petSeparateSeq))

pMeanPetSeparateDiffLowMed <- as.matrix(NA, nrow = length(petSeparateSeq))
pLoPetSeparateDiffLowMed <- as.matrix(NA, nrow = length(petSeparateSeq))
pHiPetSeparateDiffLowMed <- as.matrix(NA, nrow = length(petSeparateSeq))

pMeanPetSeparateDiffMedHigh <- as.matrix(NA, nrow = length(petSeparateSeq))
pLoPetSeparateDiffMedHigh <- as.matrix(NA, nrow = length(petSeparateSeq))
pHiPetSeparateDiffMedHigh <- as.matrix(NA, nrow = length(petSeparateSeq))

for(i in 1:length(petSeparateSeq)){
	predPetSeparateLowCong <- cbind(1,
							model2data$petUnnamedCiteLog,
							model2data$petMajCallOutLog,
							petSeparateSeq[i],
							model2data$respUnnamedCiteLog,
							model2data$respMajCallOutLog,
							model2data$respNotMajOpinLog,
							-2.909,
							model2data$pastExpertise,
							model2data$logPetNumCites,
							model2data$logRespNumCites,	
							model2data$petExperienceAdvantage,
							model2data$sgParty,
							model2data$lcDisagreement,
							model2data$amiciNet,
							model2data$helpOSG,
							model2data$netStatus,						
							model2data$oaQuestDiff,
							model2data$petUnnamedCiteLog * -2.909,
							model2data$petMajCallOutLog * -2.909,
							petSeparateSeq[i] * -2.909,
							model2data$respUnnamedCiteLog * -2.909,
							model2data$respMajCallOutLog * -2.909,
							model2data$respNotMajOpinLog * -2.909,
							issueArea,
							term,
							justice
							)
	predPetSeparateLowCong <- as.matrix(predPetSeparateLowCong)
	
	predPetSeparateMedCong <- cbind(1,
							model2data$petUnnamedCiteLog,
							model2data$petMajCallOutLog,
							petSeparateSeq[i],
							model2data$respUnnamedCiteLog,
							model2data$respMajCallOutLog,
							model2data$respNotMajOpinLog,
							0.077,
							model2data$pastExpertise,
							model2data$logPetNumCites,
							model2data$logRespNumCites,	
							model2data$petExperienceAdvantage,
							model2data$sgParty,
							model2data$lcDisagreement,
							model2data$amiciNet,
							model2data$helpOSG,
							model2data$netStatus,						
							model2data$oaQuestDiff,
							model2data$petUnnamedCiteLog * 0.077,
							model2data$petMajCallOutLog * 0.077,
							petSeparateSeq[i] * 0.077,
							model2data$respUnnamedCiteLog * 0.077,
							model2data$respMajCallOutLog * 0.077,
							model2data$respNotMajOpinLog * 0.077,
							issueArea,
							term,
							justice
							)
	predPetSeparateMedCong <- as.matrix(predPetSeparateMedCong)
	
	predPetSeparateHighCong <- cbind(1,
							model2data$petUnnamedCiteLog,
							model2data$petMajCallOutLog,
							petSeparateSeq[i],
							model2data$respUnnamedCiteLog,
							model2data$respMajCallOutLog,
							model2data$respNotMajOpinLog,
							2.918,
							model2data$pastExpertise,
							model2data$logPetNumCites,
							model2data$logRespNumCites,	
							model2data$petExperienceAdvantage,
							model2data$sgParty,
							model2data$lcDisagreement,
							model2data$amiciNet,
							model2data$helpOSG,
							model2data$netStatus,						
							model2data$oaQuestDiff,
							model2data$petUnnamedCiteLog * 2.918,
							model2data$petMajCallOutLog * 2.918,
							petSeparateSeq[i] * 2.918,
							model2data$respUnnamedCiteLog * 2.918,
							model2data$respMajCallOutLog * 2.918,
							model2data$respNotMajOpinLog * 2.918,
							issueArea,
							term,
							justice
							)
	predPetSeparateHighCong <- as.matrix(predPetSeparateHighCong)
	
	ppPetSeparateLowCong <- as.matrix(NA, row = n.draws)
	ppPetSeparateMedCong <- as.matrix(NA, row = n.draws)
	ppPetSeparateHighCong <- as.matrix(NA, row = n.draws)
	ppPetSeparateDiffLowHigh <- as.matrix(NA, row = n.draws)
	ppPetSeparateDiffLowMed <- as.matrix(NA, row = n.draws)
	ppPetSeparateDiffMedHigh <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbPetSeparateLowCong <- ilogit(predPetSeparateLowCong %*% t(sim.coef[j, ]))
		ppPetSeparateLowCong[j] <- mean(xbPetSeparateLowCong)
		xbPetSeparateMedCong <- ilogit(predPetSeparateMedCong %*% t(sim.coef[j, ]))
		ppPetSeparateMedCong[j] <- mean(xbPetSeparateMedCong)
		xbPetSeparateHighCong <- ilogit(predPetSeparateHighCong %*% t(sim.coef[j, ]))
		ppPetSeparateHighCong[j] <- mean(xbPetSeparateHighCong)
		ppPetSeparateDiffLowHigh[j] <- mean(xbPetSeparateHighCong - xbPetSeparateLowCong)
		ppPetSeparateDiffLowMed[j] <- mean(xbPetSeparateMedCong - xbPetSeparateLowCong)
		ppPetSeparateDiffMedHigh[j] <- mean(xbPetSeparateHighCong - xbPetSeparateMedCong)
		
	}
	
	pMeanPetSeparateLowCong[i] <- mean(ppPetSeparateLowCong)
	pLoPetSeparateLowCong[i] <- quantile(ppPetSeparateLowCong, 0.025, na.rm = TRUE)
	pHiPetSeparateLowCong[i] <- quantile(ppPetSeparateLowCong, 0.975, na.rm = TRUE)
	
	pMeanPetSeparateMedCong[i] <- mean(ppPetSeparateMedCong)
	pLoPetSeparateMedCong[i] <- quantile(ppPetSeparateMedCong, 0.025, na.rm = TRUE)
	pHiPetSeparateMedCong[i] <- quantile(ppPetSeparateMedCong, 0.975, na.rm = TRUE)
	
	pMeanPetSeparateHighCong[i] <- mean(ppPetSeparateHighCong)
	pLoPetSeparateHighCong[i] <- quantile(ppPetSeparateHighCong, 0.025, na.rm = TRUE)
	pHiPetSeparateHighCong[i] <- quantile(ppPetSeparateHighCong, 0.975, na.rm = TRUE)
	
	pMeanPetSeparateDiffLowHigh[i] <- mean(ppPetSeparateDiffLowHigh)
	pLoPetSeparateDiffLowHigh[i] <- quantile(ppPetSeparateDiffLowHigh, 0.025, na.rm = TRUE)
	pHiPetSeparateDiffLowHigh[i] <- quantile(ppPetSeparateDiffLowHigh, 0.975, na.rm = TRUE)
	
	pMeanPetSeparateDiffLowMed[i] <- mean(ppPetSeparateDiffLowMed)
	pLoPetSeparateDiffLowMed[i] <- quantile(ppPetSeparateDiffLowMed, 0.025, na.rm = TRUE)
	pHiPetSeparateDiffLowMed[i] <- quantile(ppPetSeparateDiffLowMed, 0.975, na.rm = TRUE)
	
	pMeanPetSeparateDiffMedHigh[i] <- mean(ppPetSeparateDiffMedHigh)
	pLoPetSeparateDiffMedHigh[i] <- quantile(ppPetSeparateDiffMedHigh, 0.025, na.rm = TRUE)
	pHiPetSeparateDiffMedHigh[i] <- quantile(ppPetSeparateDiffMedHigh, 0.975, na.rm = TRUE)
	
}

cbind(pLoPetSeparateDiffLowHigh, pMeanPetSeparateDiffLowHigh, pHiPetSeparateDiffLowHigh)
cbind(pLoPetSeparateDiffLowMed, pMeanPetSeparateDiffLowMed, pHiPetSeparateDiffLowMed)
cbind(pLoPetSeparateDiffMedHigh, pMeanPetSeparateDiffMedHigh, pHiPetSeparateDiffLowHigh)

count <- c(0, 1, 2)

petCiteSeparateXcong <- data.frame(count, petSeparateSeq, pLoPetSeparateLowCong, pMeanPetSeparateLowCong, pHiPetSeparateLowCong, pLoPetSeparateHighCong, pMeanPetSeparateHighCong, pHiPetSeparateHighCong, pLoPetSeparateMedCong, pMeanPetSeparateMedCong, pHiPetSeparateMedCong)
colnames(petCiteSeparateXcong) <- c("count", "separateCite", "lbLow", "estLow", "ubLow", "lbHigh", "estHigh", "ubHigh", "lbMed", "estMed", "ubMed")
petCiteSeparateXcong$pet <- 1

petSeparateCiteTable <- ggplot(petCiteSeparateXcong, aes(x = separateCite)) + 
				geom_line(aes(y = estLow), color = "red") +
				geom_line(aes(y = lbLow), linetype = "dashed", color = "red") +
				geom_line(aes(y = ubLow), linetype = "dashed", color = "red") +
				geom_line(aes(y = estHigh), color = "blue") +
				geom_line(aes(y = lbHigh), linetype = "dashed", color = "blue") +
				geom_line(aes(y = ubHigh), linetype = "dashed", color = "blue") +
				geom_line(aes(y = estMed), color = "green") +
				geom_line(aes(y = lbMed), linetype = "dashed", color = "green") +
				geom_line(aes(y = ubMed), linetype = "dashed", color = "green") + 
				theme_light()								
petSeparateCiteTable

#########################
### STEP 6A: RESP SEP ###
#########################

set.seed(19870302)
n.draws <- 1000
sim.coef <- coef(sim(model2logCites, n.draws))
colnames(as.data.frame(sim.coef))
sim.coef <- as.data.frame(sim.coef)
colnames(sim.coef)

respSeparateSeq <- c(0, 0.693147, 1.098612)
					
pMeanRespSeparateLowCong <- as.matrix(NA, nrow = length(respSeparateSeq))
pLoRespSeparateLowCong <- as.matrix(NA, nrow = length(respSeparateSeq))
pHiRespSeparateLowCong <- as.matrix(NA, nrow = length(respSeparateSeq))

pMeanRespSeparateMedCong <- as.matrix(NA, nrow = length(respSeparateSeq))
pLoRespSeparateMedCong <- as.matrix(NA, nrow = length(respSeparateSeq))
pHiRespSeparateMedCong <- as.matrix(NA, nrow = length(respSeparateSeq))

pMeanRespSeparateHighCong <- as.matrix(NA, nrow = length(respSeparateSeq))
pLoRespSeparateHighCong <- as.matrix(NA, nrow = length(respSeparateSeq))
pHiRespSeparateHighCong <- as.matrix(NA, nrow = length(respSeparateSeq))

pMeanRespSeparateDiffLowHigh <- as.matrix(NA, nrow = length(respSeparateSeq))
pLoRespSeparateDiffLowHigh <- as.matrix(NA, nrow = length(respSeparateSeq))
pHiRespSeparateDiffLowHigh <- as.matrix(NA, nrow = length(respSeparateSeq))

pMeanRespSeparateDiffLowMed <- as.matrix(NA, nrow = length(respSeparateSeq))
pLoRespSeparateDiffLowMed <- as.matrix(NA, nrow = length(respSeparateSeq))
pHiRespSeparateDiffLowMed <- as.matrix(NA, nrow = length(respSeparateSeq))

pMeanRespSeparateDiffMedHigh <- as.matrix(NA, nrow = length(respSeparateSeq))
pLoRespSeparateDiffMedHigh <- as.matrix(NA, nrow = length(respSeparateSeq))
pHiRespSeparateDiffMedHigh <- as.matrix(NA, nrow = length(respSeparateSeq))

for(i in 1:length(respSeparateSeq)){
	predRespSeparateLowCong <- cbind(1,
							model2data$petUnnamedCiteLog,
							model2data$petMajCallOutLog,
							model2data$petNotMajOpinLog,
							model2data$respUnnamedCiteLog,
							model2data$respMajCallOutLog,
							respSeparateSeq[i],
							-2.909,
							model2data$pastExpertise,
							model2data$logPetNumCites,
							model2data$logRespNumCites,	
							model2data$petExperienceAdvantage,
							model2data$sgParty,
							model2data$lcDisagreement,
							model2data$amiciNet,
							model2data$helpOSG,
							model2data$netStatus,						
							model2data$oaQuestDiff,
							model2data$petUnnamedCiteLog * -2.909,
							model2data$petMajCallOutLog * -2.909,
							model2data$petNotMajOpinLog * -2.909,
							model2data$respUnnamedCiteLog * -2.909,
							model2data$respMajCallOutLog * -2.909,
							respSeparateSeq[i] * -2.909,
							issueArea,
							term,
							justice
							)
	predRespSeparateLowCong <- as.matrix(predRespSeparateLowCong)
	
	predRespSeparateMedCong <- cbind(1,
							model2data$petUnnamedCiteLog,
							model2data$petMajCallOutLog,
							model2data$petNotMajOpinLog,
							model2data$respUnnamedCiteLog,
							model2data$respMajCallOutLog,
							respSeparateSeq[i],
							0.077,
							model2data$pastExpertise,
							model2data$logPetNumCites,
							model2data$logRespNumCites,	
							model2data$petExperienceAdvantage,
							model2data$sgParty,
							model2data$lcDisagreement,
							model2data$amiciNet,
							model2data$helpOSG,
							model2data$netStatus,						
							model2data$oaQuestDiff,
							model2data$petUnnamedCiteLog * 0.077,
							model2data$petMajCallOutLog * 0.077,
							model2data$petNotMajOpinLog * 0.077,
							model2data$respUnnamedCiteLog * 0.077,
							model2data$respMajCallOutLog * 0.077,
							respSeparateSeq[i] * 0.077,
							issueArea,
							term,
							justice
							)
	predRespSeparateMedCong <- as.matrix(predRespSeparateMedCong)
	
	predRespSeparateHighCong <- cbind(1,
							model2data$petUnnamedCiteLog,
							model2data$petMajCallOutLog,
							model2data$petNotMajOpinLog,
							model2data$respUnnamedCiteLog,
							model2data$respMajCallOutLog,
							respSeparateSeq[i],
							2.918,
							model2data$pastExpertise,
							model2data$logPetNumCites,
							model2data$logRespNumCites,	
							model2data$petExperienceAdvantage,
							model2data$sgParty,
							model2data$lcDisagreement,
							model2data$amiciNet,
							model2data$helpOSG,
							model2data$netStatus,						
							model2data$oaQuestDiff,
							model2data$petUnnamedCiteLog * 2.918,
							model2data$petMajCallOutLog * 2.918,
							model2data$petNotMajOpinLog * 2.918,
							model2data$respUnnamedCiteLog * 2.918,
							model2data$respMajCallOutLog * 2.918,
							respSeparateSeq[i] * 2.918,
							issueArea,
							term,
							justice
							)
	predRespSeparateHighCong <- as.matrix(predRespSeparateHighCong)
	
	ppRespSeparateLowCong <- as.matrix(NA, row = n.draws)
	ppRespSeparateMedCong <- as.matrix(NA, row = n.draws)
	ppRespSeparateHighCong <- as.matrix(NA, row = n.draws)
	ppRespSeparateDiffLowHigh <- as.matrix(NA, row = n.draws)
	ppRespSeparateDiffLowMed <- as.matrix(NA, row = n.draws)
	ppRespSeparateDiffMedHigh <- as.matrix(NA, row = n.draws)
	
	for(j in 1:n.draws){
		xbRespSeparateLowCong <- ilogit(predRespSeparateLowCong %*% t(sim.coef[j, ]))
		ppRespSeparateLowCong[j] <- mean(xbRespSeparateLowCong)
		xbRespSeparateMedCong <- ilogit(predRespSeparateMedCong %*% t(sim.coef[j, ]))
		ppRespSeparateMedCong[j] <- mean(xbRespSeparateMedCong)
		xbRespSeparateHighCong <- ilogit(predRespSeparateHighCong %*% t(sim.coef[j, ]))
		ppRespSeparateHighCong[j] <- mean(xbRespSeparateHighCong)
		ppRespSeparateDiffLowHigh[j] <- mean(xbRespSeparateHighCong - xbRespSeparateLowCong)
		ppRespSeparateDiffLowMed[j] <- mean(xbRespSeparateMedCong - xbRespSeparateLowCong)
		ppRespSeparateDiffMedHigh[j] <- mean(xbRespSeparateHighCong - xbRespSeparateMedCong)
		
	}
	
	pMeanRespSeparateLowCong[i] <- mean(ppRespSeparateLowCong)
	pLoRespSeparateLowCong[i] <- quantile(ppRespSeparateLowCong, 0.025, na.rm = TRUE)
	pHiRespSeparateLowCong[i] <- quantile(ppRespSeparateLowCong, 0.975, na.rm = TRUE)
	
	pMeanRespSeparateMedCong[i] <- mean(ppRespSeparateMedCong)
	pLoRespSeparateMedCong[i] <- quantile(ppRespSeparateMedCong, 0.025, na.rm = TRUE)
	pHiRespSeparateMedCong[i] <- quantile(ppRespSeparateMedCong, 0.975, na.rm = TRUE)
	
	pMeanRespSeparateHighCong[i] <- mean(ppRespSeparateHighCong)
	pLoRespSeparateHighCong[i] <- quantile(ppRespSeparateHighCong, 0.025, na.rm = TRUE)
	pHiRespSeparateHighCong[i] <- quantile(ppRespSeparateHighCong, 0.975, na.rm = TRUE)
	
	pMeanRespSeparateDiffLowHigh[i] <- mean(ppRespSeparateDiffLowHigh)
	pLoRespSeparateDiffLowHigh[i] <- quantile(ppRespSeparateDiffLowHigh, 0.025, na.rm = TRUE)
	pHiRespSeparateDiffLowHigh[i] <- quantile(ppRespSeparateDiffLowHigh, 0.975, na.rm = TRUE)
	
	pMeanRespSeparateDiffLowMed[i] <- mean(ppRespSeparateDiffLowMed)
	pLoRespSeparateDiffLowMed[i] <- quantile(ppRespSeparateDiffLowMed, 0.025, na.rm = TRUE)
	pHiRespSeparateDiffLowMed[i] <- quantile(ppRespSeparateDiffLowMed, 0.975, na.rm = TRUE)
	
	pMeanRespSeparateDiffMedHigh[i] <- mean(ppRespSeparateDiffMedHigh)
	pLoRespSeparateDiffMedHigh[i] <- quantile(ppRespSeparateDiffMedHigh, 0.025, na.rm = TRUE)
	pHiRespSeparateDiffMedHigh[i] <- quantile(ppRespSeparateDiffMedHigh, 0.975, na.rm = TRUE)
	
}

cbind(pLoRespSeparateDiffLowHigh, pMeanRespSeparateDiffLowHigh, pHiRespSeparateDiffLowHigh)
cbind(pLoRespSeparateDiffLowMed, pMeanRespSeparateDiffLowMed, pHiRespSeparateDiffLowMed)
cbind(pLoRespSeparateDiffMedHigh, pMeanRespSeparateDiffMedHigh, pHiRespSeparateDiffLowHigh)

count <- c(0, 1, 2)

respCiteSeparateXcong <- data.frame(count, respSeparateSeq, pLoRespSeparateLowCong, pMeanRespSeparateLowCong, pHiRespSeparateLowCong, pLoRespSeparateHighCong, pMeanRespSeparateHighCong, pHiRespSeparateHighCong, pLoRespSeparateMedCong, pMeanRespSeparateMedCong, pHiRespSeparateMedCong)
colnames(respCiteSeparateXcong) <- c("count", "separateCite", "lbLow", "estLow", "ubLow", "lbHigh", "estHigh", "ubHigh", "lbMed", "estMed", "ubMed")
respCiteSeparateXcong$pet <- 0

respSeparateCiteTable <- ggplot(respCiteSeparateXcong, aes(x = separateCite)) + 
				geom_line(aes(y = estLow), color = "red") +
				geom_line(aes(y = lbLow), linetype = "dashed", color = "red") +
				geom_line(aes(y = ubLow), linetype = "dashed", color = "red") +
				geom_line(aes(y = estHigh), color = "blue") +
				geom_line(aes(y = lbHigh), linetype = "dashed", color = "blue") +
				geom_line(aes(y = ubHigh), linetype = "dashed", color = "blue") +
				geom_line(aes(y = estMed), color = "green") +
				geom_line(aes(y = lbMed), linetype = "dashed", color = "green") +
				geom_line(aes(y = ubMed), linetype = "dashed", color = "green") + 
				theme_light()								
respSeparateCiteTable

#####################################
### STEP 6C: WRITE OUT FOR GRAPHS ###
#####################################

citesSepXcong <- rbind(petCiteSeparateXcong, respCiteSeparateXcong)

write.csv(citesSepXcong, "Analysis2SeparatePredProbs.csv")


